function glm = mdm_ComputeGLM(hfile, opts)
% MDM::ComputeGLM  - compute a GLM from an MDM file
%
% FORMAT:       glm = mdm.ComputeGLM([options])
%
% Input fields:
%
%       options     optional 1x1 struct with fields
%        .ithresh   intensity threshold, default: 100
%        .mask      optional masking, default: no mask
%        .outfile   output filename of GLM file, default: no saving
%
% Output fields:
%
%       glm         GLM object
%
% Note: if RFX flag in MDM is set to true, predictor separation will be
%       set to "Subjects".

% TODO:
%        .ar1       boolean flag, correction for AR(1), default: no

% Version:  v0.8a
% Build:    9102118
% Date:     Oct-21 2007, 6:00 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'mdm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
fdt = lower(bc.TypeOfFunctionalData(:)');

% check options
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'ar1') || ...
   ~islogical(opts.ar1) || ...
    isempty(opts.ar1)
    opts.ar1 = false;
else
    opts.ar1 = opts.ar1(1);
end
if ~isfield(opts, 'ithresh') || ...
   ~isnumeric(opts.ithresh) || ...
    numel(opts.ithresh) ~= 1
    opts.ithresh = 100;
else
    opts.ithresh = double(real(opts.ithresh(1)));
end
opts.mskf = '';
if ~isfield(opts, 'mask') || ...
    isempty(opts.mask) || ...
   (~ischar(opts.mask) && ...
    (numel(opts.mask) ~= 1 || ...
    ~isBVQXfile(opts.mask, 'msk')))
    opts.mask = [];
elseif ischar(opts.mask)
    try
        opts.mskf = opts.mask(:)';
        opts.mask = BVQXfile(opts.mskf);
        if ~isBVQXfile(opts.mask, 'msk')
            error('INVALID_MASK');
        end
    catch
        warning( ...
            'BVQXfile:BadArgument', ...
            'Invalid mask file specified: %s.', ...
            opts.mskf ...
        );
        opts.mask = [];
        opts.mskf = '';
    end
end
if opts.rfx
    opts.seppred = 2;
end

% get subject IDs
xtcfiles = regexprep(bc.XTC_RTC(:, end - 1), ...
    '^.*[\/\\]([^\/\\]*)$', '$1');
subjtokens = regexprep(xtcfiles, '^([^_]*)_.*$', '$1');
subjids = unique(subjtokens);
if any(strcmp(subjids, ''))
    error( ...
        'BVQXfile:InvalidObject', ...
        'Invalid subject IDs for some subjects.' ...
    );
end
numsubs = numel(subjids);

% determine progress bar capabilities
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', 'Computing multi-study GLM...');
    BVQXprogress(pbar, 0, 'Checking referenced files...', 'visible', 0, 1);
catch
    pbar = [];
end

% check list of files
cfs = struct('autofind', true, 'silent', true);
try
    hfile = mdm_CheckFiles(hfile, cfs);
catch
    if ~isempty(pbar)
        closebar(pbar);
    end
    error( ...
        'BVQXfile:InternalError', ...
        'Error finding referenced files.' ...
    );
end
bc = bvqxfile_getcont(hfile.L);
rfiles = bc.XTC_RTC;
rfobjs = cell(size(rfiles));
numstudy = size(rfiles, 1);

% open all files with transio
BVQXfile(0, 'transiosize', fdt, 1e5);
for fc = 1:numel(rfobjs)
    try
        rfobjs{fc} = BVQXfile(rfiles{fc});
    catch
        clearbvqxobjects(rfobjs(:));
        rethrow(lasterror);
    end
end

% check ?TCs
if ~isempty(pbar)
    BVQXprogress(pbar, 0, 'Checking dims of study 1...', 'visible', 0, numstudy);
end
numtp = 0;
try
    if fdt(1) == 'm'
        numvert = rfobjs{1, end - 1}.BVQXCONT.NrOfVertices;
        ssmvert = rfobjs{1}.BVQXCONT.NrOfTargetVertices;
        if ~isempty(opts.mask)
            opts.mask = '';
        end
        outsize = numvert;
        ststr = cell2struct(cell(0, 0, 4), ...
            {'NrOfTimePoints', 'NameOfAnalyzedFile', 'NameOfSSMFile', 'NameOfRTCFile'}, 3);
    else
        vtcres = rfobjs{1, end - 1}.BVQXCONT.Resolution;
        vtcsize = size(rfobjs{1, end - 1}.BVQXCONT.VTCData);
        vtcpos = [  rfobjs{1, end - 1}.BVQXCONT.XStart, ...
                    rfobjs{1, end - 1}.BVQXCONT.YStart, ...
                    rfobjs{1, end - 1}.BVQXCONT.XStart];
        if ~isempty(opts.mask) && ...
           (opts.mres ~= vtcres || ...
            any(opts.mpos ~= vtcpos) || ...
            any(size(opts.mask) ~= vtcsize))
            error('BAD_MASK_DIMS');
        end
        outsize = vtcsize;
        ststr = cell2struct(cell(0, 0, 3), ...
            {'NrOfTimePoints', 'NameOfAnalyzedFile', 'NameOfRTCFile'}, 3);
    end
    preds = rfobjs{1, end}.BVQXCONT.PredictorNames;
    numpred = numel(preds);
    for fc = 1:numstudy
        if ~isempty(pbar)
            BVQXprogress(pbar, fc, sprintf('Checking dims of study %d...', fc));
        end
        if fdt(1) == 'm'
            if rfobjs{fc, end - 1}.BVQXCONT.NrOfVertices ~= numvert || ...
                rfobjs{fc}.BVQXCONT.NrOfTargetVertices ~= ssmvert
                error('BAD_MTC_DIMS');
            end
            if rfobjs{fc, end - 1}.BVQXCONT.NrOfTimePoints ~= ...
                    rfobjs{fc, end}.BVQXCONT.NrOfDataPoints
                error('BAD_TIME_DIM');
            end
            ststr(fc).NameOfSSMFile = rfobjs{fc, 1}.BVQXFILE;
        else
            chksize = size(rfobjs{fc, end - 1}.BVQXCONT.VTCData);
            chkpos = [  rfobjs{fc, end - 1}.BVQXCONT.XStart, ...
                        rfobjs{fc, end - 1}.BVQXCONT.YStart, ...
                        rfobjs{fc, end - 1}.BVQXCONT.XStart];
            if any(chksize(2:4) ~= vtcsize(2:4)) || ...
                any(chkpos ~= vtcpos)
                error('BAD_VTC_DIMS');
            end
            if rfobjs{fc, end - 1}.BVQXCONT.NrOfVolumes ~= ...
                    rfobjs{fc, end}.BVQXCONT.NrOfDataPoints
                error('BAD_TIME_DIM');
            end
        end
        ststr(fc).NameOfAnalyzedFile = rfobjs{fc, end - 1}.BVQXFILE;
        ststr(fc).NameOfRTCFile = rfobjs{fc, 2}.BVQXFILE;
        spred = rfobjs{fc, end}.BVQXCONT.PredictorNames;
        if numel(spred) ~= numpred || ...
           ~all(strcmpi(spred, preds))
            error('BAD_PREDICTORS');
        end
        ststr(fc).NrOfTimePoints = rfobjs{fc, end}.BVQXCONT.NrOfDataPoints;
        numtp = numtp + ststr(fc).NrOfTimePoints;
    end
catch
    BVQXfile(0, 'transiosize' , fdt, tioosz);
    if ~isempty(pbar);
        closebar(pbar);
    end
    error( ...
        'BVQXfile:InternalError', ...
        'MDM referenced files must match in layout characteristics.' ...
    );
end

% take care of baseline
tnumpred = numpred + 1;

% create GLM object and empty contents
glm = BVQXfile('new:glm');
glm.BVQXCONT = struct;
prgb = repmat([255; 0; 0; 0], [1, 3]);

% some initial settings
glm.BVQXCONT.FileVersion = 3;
if fdt(1) == 'm'
    glm.BVQXCONT.ProjectType = 2;
else
    glm.BVQXCONT.ProjectType = 1;
end
glm.BVQXCONT.NrOfTimePoints = numtp;
snumpred = tnumpred * numsubs;
BetaMaps = single(zeros([outsize, 1]));
if opts.rfx
    glm.BVQXCONT.ProjectTypeRFX = 1;
    glm.BVQXCONT.NrOfSubjects = numsubs;
    glm.BVQXCONT.NrOfSubjectPredictors = tnumpred;
    RFXGlobalMap = single(zeros([outsize, 1]));
    rnumpred = tnumpred;
    BetaMaps(:, :, :, 1:rnumpred, 1:numsubs) = 0;
else
    glm.BVQXCONT.ProjectTypeRFX = 0;
    switch (opts.seppred)
        case {0}
            rnumpred = tnumpred;
        case {1}
            rnumpred = tnumpred * numstudy;
        case {2}
            rnumpred = tnumpred * numsubs;
    end
    snumpred = rnumpred;
    MultipleRegressionR = single(zeros([outsize, 1]));
    MCorrSS = single(zeros([outsize, 1]));
    BetaMaps(:, :, :, 1:tnumpres) = 0;
    XY = single(zeros([outsize, 1]));
    XY(:, :, :, 1:tnumpres) = 0;
    TimeCourseMean = single(zeros([outsize, 1]));
end
glm.BVQXCONT.NrOfPredictors = snumpred;
glm.BVQXCONT.NrOfStudies = numstudy;
glm.BVQXCONT.SeparatePredictors = opts.seppred;
glm.BVQXCONT.TransformationType = opts.norm;
glm.BVQXCONT.SerialCorrelation = opts.ar1;
glm.BVQXCONT.MeanAR1Pre = 0;
glm.BVQXCONT.MeanAR1Post = 0;
if fdt(1) == 'm'
    glm.BVQXCONT.Resolution = 3;
    glm.BVQXCONT.NrOfVertices = numvert;
else
    glm.BVQXCONT.Resolution = rfobjs{1, end - 1}.BVQXCONT.Resolution;
    glm.BVQXCONT.XStart = chkpos(1);
    glm.BVQXCONT.XEnd   = rfobjs{1, end - 1}.BVQXCONT.XEnd;
    glm.BVQXCONT.YStart = chkpos(2);
    glm.BVQXCONT.YEnd   = rfobjs{1, end - 1}.BVQXCONT.YEnd;
    glm.BVQXCONT.ZStart = chkpos(3);
    glm.BVQXCONT.ZEnd   = rfobjs{1, end - 1}.BVQXCONT.ZEnd;
end
if isempty(opts.mask)
    glm.BVQXCONT.CortexBasedStatistics = 0;
    glm.BVQXCONT.CortexBasedStatisticsMaskFile = '';
    opts.mask = true([outsize, 1]);
else
    glm.BVQXCONT.CortexBasedStatistics = 1;
    glm.BVQXCONT.CortexBasedStatisticsMaskFile = opts.mskf;
end
glm.BVQXCONT.Study = ststr;
glm.BVQXCONT.Predictor = struct;

% generating FFX design matrix
prstr = cell2struct(cell(0, 0, 3), {'Name1', 'Name2', 'RGB'}, 3);
fpnum = 1;
rpr = rfobjs{1, end}.BVQXCONT.Predictornames;
if ~opts.rfx 
    if ~isempty(pbar)
        BVQXprogress(pbar, 0, 'Creating FFX design matrix...', 'visible', 0, 1);
    end
    XX = zeros(numtp, rnumpred);
    ctp = 1;
    cpr = 1;
    switch (opts.seppred)
        case {0}
            for spc = 1:numel(rpr)
                prstr(fpnum).Name1 = sprintf('Predictor: %d', fpnum);
                prstr(fpnum).Name2 = rpr{spc};
                prstr(fpnum).RGB = prgb;
                fpnum = fpnum + 1;
            end
            for sc = 1:numstudy
                rtc = rfobjs{sc, end}.BVQXCONT.RTCMatrix;
                rsz = size(rtc);
                XX(ctp:(ctp + rsz(1) - 1), 1:rsz(2)) = rtc;
                XX(ctp:(ctp + rsz(1) - 1), (rsz(2) + sc)) = 1;
                ctp = ctp + rsz(1);
                prstr(fpnum).Name1 = sprintf('Predictor: %d', fpnum);
                prstr(fpnum).Name2 = 'Study signal level (confound)';
                prstr(fpnum).RGB = prgb;
                fpnum = fpnum + 1;
            end
        case {1}
            for sc = 1:numstudy
                rtc = rfobjs{sc, end}.BVQXCONT;
                rpr = rtc.PredictorNames;
                rtc = rtc.RTCMatrix;
                rsz = size(rtc);
                XX(ctp:(ctp + rsz(1) - 1), cpr:(cpr + rsz(2) - 1)) = rtc;
                XX(ctp:(ctp + rsz(1) - 1), (end - numstudy + sc)) = 1;
                ctp = ctp + rsz(1);
                cpr = cpr + rsz(2);
                for spc = 1:numel(rpr)
                    prstr(fpnum).Name1 = sprintf('Predictor: %d', fpnum);
                    prstr(fpnum).Name2 = sprintf('Study %d: %s', sc, rpr{spc});
                    prstr(fpnum).RGB = prgb;
                    fpnum = fpnum + 1;
                end
            end
            for sc = 1:numstudy
                prstr(fpnum).Name1 = sprintf('Predictor: %d', fpnum);
                prstr(fpnum).Name2 = 'Study signal level (confound)';
                prstr(fpnum).RGB = prgb;
                fpnum = fpnum + 1;
            end
        case {2}
            subjtoken = subjtokens{1};
            for sc = 1:numstudy
                rtc = rfobjs{sc, end}.BVQXCONT;
                rpr = rtc.PredictorNames;
                rtc = rtc.RTCMatrix;
            end
    end
    if ~isempty(pbar)
        BVQXprogress(pbar, 0.25, 'Inverting FFX design matrix...');
    end
    iXX = pinv(XX);
    glm.BVQXCONT.Predictor = prstr;
end

% calc depends of Separate predictors / RFX
switch (opts.seppred + 10 * opts.rfx)
    
    % no separation (full FFX)
    case {0}
        
    % study separation (FFX with full separation)
    case {1}
        
    % subject separation (RFX-like separation but FFX model)
    case {2}
        
    % RFX
    case {12}
        
        % iterate over subjects
        for sc = 1:numsubs
        end
        
end

% fill GLM
if opts.rfx
    glm.BVQXCONT.GLMData.RFXGlobalMap = RFXGlobalMap .* opts.mask;
    glm.BVQXCONT.GLMData.BetaMaps = BetaMaps;
else
    glm.BVQXCONT.GLMData.MultipleRegressionR = ...
        MultipleRegressionR .* opts.mask;
    glm.BVQXCONT.GLMData.MCorrSS = MCorrSS .* opts.mask;
    glm.BVQXCONT.GLMData.BetaMaps = BetaMaps;
    glm.BVQXCONT.GLMData.XY = XY;
    glm.BVQXCONT.GLMData.TimeCourseMean = TimeCourseMean .* opts.mask;
end

% re-set transio setting
BVQXfile(0, 'transiosize', fdt, tioosz);

% close pbar
if ~isempty(pbar)
    closebar(pbar);
end
